#include "MITC3PlateCalculator.h"

#include<iostream>
using SG::DataStructure::Common::Real;
using SG::Algebra::Matrix;
using SG::FEMSolver::Element::ShellSectionProp1;


std::size_t MITC3PlateCalculator::s_nNode = 3;

void MITC3PlateCalculator::Initialize ()
{
}

/// @brief 应变矩阵计算
void MITC3PlateCalculator::ComputeB (std::size_t iGauss, const Matrix& coord,  Matrix& _OUT Bs, _OUT Matrix& Bb, _OUT Real& j) const
{
    // 弯曲部分使用长应变
    //行列式求面积
    Real twoA = coord(1, 0) * coord(2, 1) + coord(2, 0) * coord(0, 1) + coord(0, 0) * coord(1, 1) -
                    coord(1, 0) * coord(0, 1) - coord(2, 0) * coord(1, 1) - coord(0, 0) * coord(2, 1);

    j = twoA * 0.5;

    Bb = Matrix(3, 9);
    // node1
    const Real dN1dx = - (coord(1, 1) - coord(2, 1)) / twoA;
    const Real dN1dy = - (coord(2, 0) - coord(1, 0)) / twoA;
                     Bb(0, 2) = -dN1dx; 
    Bb(1, 1) = dN1dy;
    Bb(2, 1) = dN1dx; Bb(2, 2) = -dN1dy;
    // node2
    const Real dN2dx = - (coord(2, 1) - coord(0, 1)) / twoA;
    const Real dN2dy = - (coord(0, 0) - coord(2, 0)) / twoA;
                     Bb(0, 5) = -dN2dx;
    Bb(1, 4) = dN2dy;
    Bb(2, 4) = dN2dx; Bb(2, 5) = -dN2dy;
    // node3
    const Real dN3dx = - (coord(0, 1) - coord(1, 1)) / twoA;
    const Real dN3dy = - (coord(1, 0) - coord(0, 0)) / twoA;
                     Bb(0, 8) = -dN3dx;
    Bb(1, 7) = dN3dy;
    Bb(2, 7) = dN3dx; Bb(2, 8) = -dN3dy;

    /** @brief 使用 Mix Interpolation of Tensor Component技术克服剪切锁定现象
     * 
     *  使用边中点处切应力对积分点切应变插值 [既增强假定应变（Enhanced Assumed Strain）技术]
     */
    const auto xi = m_integForm.GetPoint(iGauss).m_x;
    const auto eta = m_integForm.GetPoint(iGauss).m_y;
    
    Bs = Matrix(2, 9);
    Real x12 = coord(1, 0) - coord(0, 0); Real y12 = coord(1, 1) - coord(0, 1);

    Real x31 = coord(0, 0) - coord(2, 0);
    Real y31 = coord(0, 1) - coord(2, 1);
    Real x13 = - x31;
    Real y13 = - y31; 

    Real x23 = coord(2, 0) - coord(1, 0);
    Real y23 = coord(2, 1) - coord(1, 1);

    //node 1
    Bs(0, 0) = -2.0;                Bs(1, 0) = -2.0; 
    Bs(0, 1) = -x12 - x23 * eta; Bs(0, 2) = -y12 - y23 * eta;
    Bs(1, 1) = -x13 + x23 * xi;  Bs(1, 2) = -y13 + y23 * xi;
     
    //node 2
    Bs(0, 3) = 2.0;                 Bs(1, 3) = 0.0;
    Bs(0, 4) = -x12 -x31 * eta; Bs(0, 5) = -y12 - y31 * eta;
    Bs(1, 4) =       x31 * xi;  Bs(1, 5) =        y31 * xi;
     
    //node 3
    Bs(0, 6) = 0.0;                 Bs(1, 6) = 2;
    Bs(0, 7) =      - x12 *eta;     Bs(0, 8) =      - y12 * eta; 
    Bs(1, 7) = -x13 + x12 * xi;     Bs(1, 8) = -y13 + y12 * xi;
            
 
 
    Matrix J = Matrix(2, 2, {x12, y12, x13, y13});   
    Bs = inverse2X2(J) * (Bs * 0.5);
 
    Matrix T = Matrix(9, 9);
    T(0, 0) = 1; T(1, 2) = -1; T(2, 1) = 1;
    T(3, 3) = 1; T(4, 5) = -1; T(5, 4) = 1;
    T(6, 6) = 1; T(7, 8) = -1; T(8, 7) = 1;
    
    Bs = Bs * T;// 初等变换交换列

}


// void MITC3PlateCalculator::ComputeBs (std::size_t iGauss, Real j, const Matrix& coord,  Matrix& _OUT Bs) const
// {

// }

Matrix MITC3PlateCalculator::ComputeStiffness (const Matrix& coord, const Real* const thickness, const ShellSectionProp1& secConstant) const
{
    Matrix ke (9, 9);
    // SecConstant.GetCb (thickness[0], _OUT Cb);
    const auto& Cb6 = secConstant.m_D;
    Matrix Cb (3,3, {
        Cb6(0,0), Cb6(1,0), Cb6(2,0),
        Cb6(1,0), Cb6(3,0), Cb6(4,0),
        Cb6(2,0), Cb6(4,0), Cb6(5,0)});
    auto IOvr12 = secConstant.FI2 /12.0;
    Cb = Cb * (thickness[0] * thickness[0]* thickness[0] * IOvr12);

    // 弯曲部分
    Matrix Bb;
    Matrix Bs;
    Real detj;
    
    // 仅仅Bb 有效
    ComputeB (0 /* 使用单点积分 */, coord, Bs, _OUT Bb, _OUT detj);
    ke += Bb.Transpose() * Cb * Bb * detj;
    
    Matrix Cs;
    const auto nGauss = m_integForm.GetPointSize();
    for (std::size_t iGauss = 0; iGauss < nGauss; ++iGauss)
    {
        // SecConstant.GetCs (thickness[iGauss], _OUT Cs);
        auto factor = secConstant.m_ts * thickness[iGauss];
        Cs = Matrix (2, 2, {
            secConstant.m_Cs(0, 0) * factor, secConstant.m_Cs(1, 0) * factor,
            secConstant.m_Cs(1, 0) * factor, secConstant.m_Cs(2, 0) * factor
        });

        // ComputeBs (iGauss, detj, coord, _OUT Bs);
        // 仅仅使用Bs
        ComputeB (iGauss, coord, Bs, _OUT Bb, _OUT detj);
        const Real tmp = detj * m_integForm.GetPoint(iGauss).m_weight; 
        ke += Bs.Transpose() * Cs * Bs * tmp;
    }
    return ke;
}